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Abstact: We studied numerically the validity of the fluctuation theorem introduced 
in [1][2] for a 2-dimensional system of particles maintained in a steady shear flow 
by Maxwell daemon boundary conditions [Sj. The theorem was found to hold if one 
considers the total phase space contraction a occuring at collisions with both walls: 
a = a 1 * +o~^ . An attempt to extend it to more local quantities and , corresponding 
to the collisions with the top or bottom wall only, gave negative results. The time decay 
of the correlations in crTd was very slow compared to that of a. 

Keywords: dynamical system, time reversibility, fluctuation theorem, local phase 
space contraction. 



1. Introduction 

The microscopic structure of systems through which there is transport of energy or momentum, 
is a central problem in non-equilibrium statistical mechanics. For gases a partial answer to this 
question (on the mesoscopic/kinetic level) is provided by the stationary solution of the Boltzmann 
equation with suitable, e.g. Maxwellian, boundary conditions. Going beyond kinetic theory has 
proven to be very difficult and we still lack a full understanding of stationary non-equilibrium 
states (SNS) on the microscopic level. 

This gap goes beyond that of computational complexity or even of technical difficulty in proving 
the validity of the formulas derived formally, something already present also in the kinetic theory 
and in equilibrium statistical mechanics. What is still missing, at the present time, are well defined 
formal procedure which would, at least in principle, provide answers to questions of physical 
relevance for SNS of macroscopic systems. Thus, we have no "statistical mechanical formula" for 
the decay of spatial or time displaced correlation functions in a SNS. In particular, we have no a 
priori formalism for deciding on the slow power law decay of the spatial correlations predicted by 
computer simulations, fluctuating hydrodynamics and confirmed by experiments [4]. We also have 
no formula for computing the stationary heat flux through a metal or plasma or the momentum flux 
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through a fluid which goes beyond linear response theory. Such a formula should, like the Green- 
Kubo formula for linear transport, depend only on on the internal Hamiltonian of the system and 
the impressed macroscopic constraints driving the system, e.g. a specified temperature or velocity 
field on the boundaries of the system. 

While the usual modeling of such systems has been via stochastic boundaries [5] there has been 
much interest recently in the study of SNS of particle systems evolving via entirely deterministic 
dynamics in the hope that dynamical systems theory will provide new insights into non-equilibrium 
behavior [6] [7] . As is well known the existence of such SNS in finite systems is incompatible with 
the usual Hamiltonian evolution, believed to accurately describe the dynamics of systems in which 
quantum effects are unimportant. For such dynamics the only realizable stationary states are those 
which depend solely on the global constants of motion. In realistic systems these are just the total 
energy and (under suitable boundary conditions) the total momentum and angular momentum. 

Adding external (globally non-gradient) forces to the dynamics, e.g. a uniform field in a system 
with periodic boundary conditions (such as would arise from a changing magnetic flux crossing 
the surface bounded by a conductor) will typically result in the system gaining energy continu- 
ously from the field, excluding the possibility of SNS. This necessitates the use of non phase space 
volume conserving forces for modeling SNS. Various models of such dynamics have been inves- 
tigated through computer simulations and heuristic analysis [7], and there have also been some 
mathematical results motivated by physical considerations of such SNS [2]. There is however, as 
already mentioned, still a very wide gap between the mathematical and physical results obtained 
from these models and a full theory of macroscopic SNS. 

While it is not clear at present whether and how the dynamical system approach will answer 
such questions it seems reasonable to explore the behavior of such SNS. This is particularly so for 
models in which the imposed, model dependent, driving and thermostating terms are confined to 
the boundaries of the system while the dynamics in the interior of the system remain realistically 
Hamiltonian. Several such models of shear flow were investigated via computer simulations and 
some heuristic analysis in [3]. Here we continue our investigation focusing on the behavior of 
the fluctuation in the phase space volume contraction a which occurs at the boundaries, a was 
found in [3] to be approximately equal, when the size of the system is sufficiently large for it to 
be in local thermal equilibrium, to the hydrodynamic entropy production inside the system. The 
same quantity, in a different model of shear flow, was studied in [1] where an interesting relation 
for its large deviations was experimentally found. This relation is now a rigorous result for large 
deviations of phase space contraction under suitable conditions on the dynamics [2] . Whether the 
results hold when the conditions are not satisfied exactly is a question of great relevance. The 
exploration of this and related quantities is the subject of the present work. 

In the next section we describe the model and the check of the Gallavotti-Cohen result [2] for 
our system. In section 3 we investigate fluctuations and time dependent correlations of the volume 
contractions and produced by collisions with the top and bottom walls. 
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2. Numerical simulations: the fluctuation theorem 

A: The model 

We consider the 2-dimensional shear flow model first introduced in [3] and further studied in [8]: 
iV identical disks of radius r evolve in the interior of the system according to Hamiltonian dynamics 
with hard core interaction among them while Maxwell deamons at the walls drive the system away 
from equilibrium. More specifically the particles move on the surface of a cylinder (i.e. the system 
has periodic boundary conditions on the vertical sides) with reflection rules at the top and bottom 
walls simulating macroscopic moving walls. The rules are as follows: when a particle collides with 
the upper wall, making an angle between the positive x-direction and the incoming velocity then 
the outgoing velocity angle tp will be given by a (to be specified) reflection rule /: 

V = /(0) (2.1) 

A similar rule applies on the lower wall with the only difference that ip and <p refer to the angle 
between the incoming and outgoing velocity and the negative x-direction. Since the modulus of 
the velocity is preserved during a collision the total kinetic energy of the system is a constant of 
the motion. 

We will assume that / is "time reversible", i.e. it satisfies = f(n — f{jr — </>)). Observe that if 
this condition holds the system is time reversible: if one inverts the velocities of all the particle 
the system traces back its past trajectory. In the numerical simulations we will consider one of the 
reflection rules introduced in [3], namely: 

V> = (tt + b) - V(7r + b) 2 - <P(<j) + 2b) (2.2) 

where b is used to modulate the intensity of the shear with b = oo representing elastic collisions 
(see Fig. 1). 




Fig. 1: Schematic representation of the dynamics of the system 

During a collision with the walls the Liouville measure is not conserved. In fact at each "reflection 
there is a phase space contraction a equal to [3]: 
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Q = _l 0g (^/W) (2.3) 
\ sin <p J 

It is natural to consider the collisions with the walls as timing events for the system, i.e. we 
can consider the function T that associates to any point X in the phase space of the system, 
i.e. the energy surface E, the first time at which a particle collides with a wall and define the 
map, S : E — > E, S(X) = F($(X,T(X))), E is the set of points at which at least one particle is 
colliding with the walls where t) is the point X t obtained from X via the flow generated by 
the dynamics and F : E — > E is the collision map, i.e. the map that change the velocity of the 
colliding particle according to the reflection rule /. 

Given a point X we define 

[r/2]-l 

a T (X)= £ a(S*(X)) (2.4) 

i=-[r/2] 

and 

M*) = ^ (2-5) 

where [t] is the biggest integer less that t, (•) + refer to the mean with respect to the forward 
invariant distribution. Call 7r T (p) the distribution function of cr T , and let 

x T (p) = T ^log^- (2.6) 
(a T )+ 7r T (-p) 

then the chaotic hypothesis of [2] implies that 

lim x T (p) = p (2-7) 



To check numerically this relation numerical simulations were carried out on systems of N = 
10, 20, 30 and 40 particles of radius R = 1 in a "square" box of size L = y/dN, where d = N/L 2 is 
the number density kept constant to d = 3.4 x 10 -2 , and with the parameter b = 25. The initial con- 
figuration was chosen randomly in such a way that the energy per particle eo = ^ ^iLo v i = 0-5. 
In the next table we give some of the interesting dynamical quantities associated with the system. 





10 


20 


30 


40 


d 


3.4 x 10" 2 


3.4 x 10" 2 


3.4 x 10" 2 




to 


20.6 


30.6 


38.1 


44.4 


V 


2.11 


3.19 


4.00 


4.67 


(a) + 


1.24 x 10" 2 


1.68 x 10" 2 


2.00 x 10" 2 


2.26 x 10" 2 


^max 


0.92 


0.74 


0.63 


0.57 



Table 1: static and dynamic quantity for the simulation where p.s.c stands for phase space contraction. 
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Here to is the mean time between successive collisions with the walls for a given particle, v is the 
number of binary collisions (i.e. collisions between two particles) between two consecutive collisions 
with the walls for a given particle and (a)+ in the mean phase space contraction rate. Finally we 
give the maximum Lyapunov exponent of the map S. The data, given without error estimates, are 
intended to give a rough image of the dynamics. 

B: The fluctuation law 

The check of the fluctuation theorem was done as in [9]. A long trajectory of 5 x 10 8 collisions 
with the walls was simulated and the phase space contraction was recorded for every 100 collisions. 
The main difference with [9] is that we did not attempt to decorrelate the adjacent data segment 
by leaving out a fixed number of collision between them. In fact in [9] this was possible because 
the self-correlation of the phase space contraction was decaying rapidly enough to leave out just 
few collisions. In the present system we find that while the total phase space contraction rate 
has rapidly decaying correlation the partial ones (to be defined precisely and studied in the next 
section) have a very slow decay. Hence we would have to discard too many collisions to decorrelate 
the adjacent data segments and we therefore decided to discard no collision also for the total phase 
space contraction to have a consistent analysis. 

In Fig. 2 we show the graph of x T (p) for N = 20 and several value of r. 

As we can see already for r = 400 the agreement between the theoretical prediction and the 
experiment is very good. 

To better follow the behavior of the fluctuation we have constructed the function x T {p) for 
N = 10, 20, 30 and 40 for several value of r from 1 to 1000. We have then used a least square fit 
to fit the experimental data with a law of the form x T (p) = XtP- We chose a one parameter fit 
because x T (0) = 0. The result are shown in Fig. 3. As can be easily seen the evaluated \ T contain 
1 within their error-bars starting from quite small value of r. Moreover, one appears to be the 
asympotical value of \ T . 

The analysis of the fluctuation law requires the construction of the whole distribution function 
7r T (p). But as r grows big fluctuations become more and more improbable and it is impossible to 
construct x T (p) for r > 1000. To go further with r we can observe that, starting with r ~ 5, the 
distribution n T (p) look very much like a Gaussian. Observe that if we assume that the system is 
Anosov the central limit theorem implies a Gaussian behavior near the maximum of the distribution 
but the observed agreement goes, in our opinion, beyond the prediction of the central limit theorem. 
Nevertheless it must be noted that the distribution can't be Gaussian because it is easy to observe 
that |a T (X)| is bounded. Fig 4 show the comparison with a Gaussian for N = 20 and several 
values of r. 

Assuming that tt t is a Gaussian immediately implies that x T (p) is linear in p and, calling C T the 
covariance of a T , i.e. C T = (a^) — l,we have that 
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Fig. 2: Graphs of x T (p) for r = 100,400,600,800 for iV = 20. The arrows represents the value 1 plus or minus the 
standard deviation of o>. 



The covariance C T can be easily computed from the data and permit us to go beyond the limit 
of t = 100 met before. As a further check we show the value of Xt as computed from the best fit 
and from the Gaussian hypothesis for small value of r in Fig. 5. 

The agreement is very good and justifies to use of eq.(2.8) for large value of r. The evaluated 
behavior of Xt for large value of r from the Gaussian hypothesis is shown in Fig. 6. 

In all cases the values of Xt are very close to 1 and we can therefore say that the prediction of 
the fluctuation law appears to be verified by our numerical simulations. 

Finally observe that the approach to 1 of Xt can be connected to the decay property of the 
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Fig. 3: Behavior of Xt as a function of r for AT = 10, 20, 30 and 40. 



auto-correlation D(t) = (a(S t (-))a(-)} — 1. In fact we have: 

C(r) = H £ - A ^ (2.9) 

t= — T t= — T 

The fast approach of % T to its limit can thus be interpreted as a rapid decay of the correlation 
for the phase space contraction of the system. We will see that this behavior changes greatly when 
we consider partial phase space contractions. 



3. Local fluctuations 

An interesting question is whether one can give a local version of the fluctuation theorem. To 
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Fig. 4: Comparison between the distribution n T (p), for N = 20 and several values of r, and a Gaussian with the same 
variance and mean. The error-bars are smaller than the dimension of the points. 

this extent observe that we can write a T (X) = a\{X) + a\.{X) where 

r/2 

al(X)= a(5 4 (X))x T (5«(X)) (3.1) 

i=-r/2 

and x ( x ) = 1 if -X" correspond to a collision of a particle with the upper wall and otherwise. An 
analogous definition holds for a\.(X). 
The generalized version of the fluctuation theorem of [10] gives 
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lim-i-log ff<*'*> = P1+P2 (3.2) 
— (al)+ 7il l (-p u -p 2 ) 

where 7r^(pi,p2) is the joint distribution of crl(X) = ot\(X) / [et\) + and cr^(X) = a\.{X)/ {a\) + . 
Due to the symmetry of the problem the two variables a\(X) and a\.(X) can be assumed to be 
identically distributed 1 . It is easy to see that if we suppose them independent, both of them will 
separately satisfy the fluctuation theorem: 

lim x\(p) = p (3-3) 

where 

t 1 ttHv) 

xl(p) = log -pV>- (3.4) 



Although this is well verified numerically it is not evident form a theoretical point of view. It is in fact easy to see that 
if the particles do not interact, i.e. they do not collide, the mean momentum of the center of mass is, generically, not zero 
and this will probably create asymmetry between the two walls. This phenomenon is probably destroyed by the "mixing" 
behavior generated by the binary collisions. It must be nevertheless noted that for small systems, like ours, the fluctuations 
of the center of mass momentum are quite big and its correlations decay slowly 
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Fig. 6: Value of %r obtained from the Gaussian hypothesis for large r 



This can be checked as we did for a T (x). In Fig. 7 we show the graph of x\{p) for N = 20 and 
several values of r. 

Observe that here we are able to construct the distribution function n T up to r = 2000. This is 
not surprising considering that cv T {X) is, in the mean, the sum of r/2 non zero terms so that it 
can be expected to have fluctuation roughly similar to those of a T / 2 {X). 

As before we can look at the best fit of x\{jp) of the form x\{p) — xIp- The results are shown in 
Fig. 8. 

It is clear that for N = 10 and 20 the fit goes significatively below 1 at a value of r for which \ T 
is very near its theoretical value of 1. Moreover xX does not seem to have reached a limiting value 
while Xt reaches its limiting value of one quite soon. 

The cases N = 30 and 40 are less clear because we are unable to construct the distribution 
function for r large enough to have a clear idea of the limiting value of Xt- But it seems reasonable 
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Fig. 7: Graphs of x' T {p) for r = 100,800,1400,2000 for TV = 20. The arrows represents the value 1 plus o minus the 
standard deviation of uj. 

to deduce from the graph that also in this case \ T will became smaller that 1. 

Considering the argument at the beginning of this section it is interesting to look at the cross 
correlation between a\ and a\. given by: 

Cl l = - 1 

whose graph is plotted in Fig. 9. 

It is interesting to note that, if we disregard the value of for small value of r, we have behavior 
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Fig. 8: Behavior of xl as a function of r for TV = 10, 20, 30 and 40. 



of is linked to that of Xt- In particular we can see that when = we have % r = 1. This 
is particularly clear in the case N = 10 and 20. 

This suggests that the higher order correlations of and are small compared to the self- 
correlation. Moreover the distribution function n\ is again very well approximated by a Gaussian 
as can be seen in Fig. 10. 

As in the previous section we can assume that the pair (a\a^) represents a bidimensial Gaussian 
variable with mean (1,1) and covariance matrix: 

Cr = (cjt %) (3 - 5) 

where 
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Fig. 9: Behavior of C± l as a function of r for N = 10, 20, 30 and 40. 



CV T = <(^) 2 >-1 C r i = <(^) 2 >-1 (3-6) 

Form now on we will suppose that Cj = C\. This is well verified in the experiment and seems 
natural from the symmetry properties of the system (see however footnote 1). With the above 
Gaussian hypothesis we can again try to follow the behavior of xl f° r large value of r using the 
expression: 

xl 

Fig. 11 shows the results of this evaluation. 
Finally we observe that we can write: 

Cl l = 2C T - C\ (3.8) 
Hence, if we assume gaussianity and take r large enough that we can consider \ T = 1, gives 
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Fig. 10: comparison between the distribution 7rj!(p) and a Gaussian with the same variance and mean. The errorbar are 
smaller than the dimension of the points. 



K T = xl-l = -= r (3-9) 

We can conclude this discussion by observing that, for finite N, the fluctuation law seems to be 
invalid if we take in consideration only half of the entropy production. Nevertheless it is interesting 
to note that x\.'^{p) still look linear and that the slope x T seems to reach a finite limit as r goes 
to infinity. Moreover this limit is very small for small iV and seems to increase with N . It would 
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Fig. 11: Behavior of xl as a function of r for N = 10, 20, 30,40. 



be intersting to fit the curve in Fig. 11 to give an estimate of the limiting value of xl as T ~^ 00 • 
We observe that, as for Xt, we can relate the approach to a limiting value of xl with the decay 
properties of (t) = (cr^(5 ,t (-))cr^(-)) — 1. In fact an equation identical to (2.9) holds for C\ with 
D^(t) in the place of D(t). A comparison of Fig. 6 and Fig. 11 immediately shows that the decay 
of correlations of 0^ is much slower that that of a. If we assume that £> T (t) ~ t~P, with (3 > 0, we 
get that: 

+ / It 

T = T . , * 

Xt Xoo "t ^ -r T /3-i 

If f3 is smaller than 2 the third term dominate, in the asymptotic behavior on the second one. 
One can try to fit (3 looking at the log-log plot of xl an d then using a least square fit to find the 
other coefficient. Unfortunately the log-log plot does not give a precise answer, so we can only say 
that (3 is probably less than 2. This implies that the approach to the limit of xl is ver Y slow and 
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its limiting value change considerably depending which (3 one uses. At the end we are unable to 
give a quantitative estimate of xb} an d we can just say that it is non zero and increases with N. 

4. Conclusion 

Our numerical results show that the Fluctuation Theorem, eq.(2.5), is well verified if we consider 
the total phase space contraction a of the model described in sec. 1 sub A. For r small enough 
we were able to construct the distribution ir T (p) and check directly the validity of eq. (2.5). For 
higher r we used a Gaussian hypothesis to compute the slope Xr of x T (p) obtaining again a very 
good agrrement with the theoretical prediction. We observe that if we interpret the Gaussianity of 
7r r (p) as a Central Limit Theorem effect we get that x T (p) is linear but no conclusion con be drawn 
on the slope % r . The validity of the fluctuation theorem predictions, toghether with the gaussian 
hypothesis, imply a Green-Kubo formula out of equilibrium as can be seen by combining eq. (2.8) 
and eq. (2.9) (see also [9]). In this sense we can think that our results depend in part on the fact 
that we are not very far from equilibrium. Increasing the shear make the experiment much harder 
becouse the probability of observing negative fluctuations decreases very rapidly with the shear. 
Moreover we expect that the gaussianity of n T (p) will be destroyed. We think, however, that if the 
shear is not too big, i.e. if the attractor can still be considered dense in phase space (see [11]), then 
the fluctuation theorem will holds (see [12] for a situation in wich 7r T (p) is clearly not Gaussian). 

We also checked the validity of the fluctuation theorem for the partial phase space contraction 
a' and crK In this situation we see that is still well aproximated by a Gaussian so that x\(p) 
appears to be linear. The slope xX however behave in a different way from what one would expect. 
Two interesting features of this behaviour are: 

1) It saturates very slowly. The fact that we considered, in the definition of a fixed number 
of collisions r with both walls clearly create a negative correlation between a\ and a\. This 
correlation can be roughly estimated as r _1 so that it probably cannot account for the slow 
behavior found in sec. 3. That behavior shows a strong correlation between the two wall 
probably due to the slow decay of fluctuation the density near the walls. We do not have a clear 
understanding of this phenomenon and we hope to come back on it in a forthcaming work. 

2) The limiting value xlo * n n °t easy to estimate but is clearly smaller than 1. Although it seems 
to increase with the number of particle we do not have enough data to draw any conclusion 
on its behaviour with N. Nevertheless the existence of a limiting non-zero value suggests the 
validity of a "local fluctuation theorem" in wich the slope 1 is replaced by some other value. This 
behavior may have relation to the recent experimental result on the local entropy production 
(supposed to be the equivalent of a phase space contraction) in a fluid moving in a convective 
cell [13]. They find there a clear linear behaviour but the slope seems to be different from one. 
The rusults are, however, still unclear and no real comparison can be made. 
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